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Q Abstract 
We consider the late time behaviour of non-tilted perfect fluid Bianchi VIIo models when the source 
l/") ' is a radiation fluid, thereby completing the analysis of the Bianchi VIIo models initiated by Wainwright 

et al in a recent paper. The models exhibit the phenomena of asymptotic self-similarity breaking and 
Weyl-curvature dominance at late times. The late time dynamics of the VIIo perfect fluid models, and in 
particular that of the radiation fluid, is a prime example of the complexity inherent in the field equations 
of general relativity. 
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1 Introduction 

on: 

ON ■ The late time behaviour of Bianchi type VIIo models with a non-tilted perfect fluid source has been shown 

' to exhibit two new dynamic features at late times, namely, the breaking of asymptotic self-similarity and 

q-i Weyl curvature dominance (see Wainwright et al (1999)[], pages 2587-8). These models therefore serve as 

counterexamples to the notion that any Bianchi universe is approximated by a self similar Bianchi universe 
at late times. The breaking of asymptotic self-similarity is characterized by oscillations in the dimensionlcss 
shear scalar that become increasingly rapid in terms of the cosmological clock time t as t — > oo. This 
behaviour also leads to the Weyl curvature dominance, which refers to the fact that certain expansion- 
' normalized scalars formed from the Weyl tensor become unbounded as r — * +oo . 

It is customary in cosmology to describe the source of the gravitational field in terms of a perfect fluid, 
for which the energy-momentum tensor has the form 

T ab = flU a U b + p (g a b + U a Ub) , (1) 

where /i > is the fluid energy density, p > the fluid pressure, and u a the fluid 4-velocity. In WHU, the 
linear equation of state 

P = (7-l)/*, (2) 

where the constant 7 satisfies | < 7 < 2, was studied. The two cases of primary interest in cosmology are 
dust (7 = 1), describing the matter-dominated epoch of the universe, and radiation (7 = ~), describing 
the radiation-dominated epoch of the early universe. The late-time behaviour^] of a radiation model would 
describe the final stages of this epoch, prior to the time at which matter and radiation decouple. In WHU 
it was shown that for Bianchi VIIo models a physically important bifurcation occurs at the value 7 = | , 
namely that the dimensionless shear scalar £ (see Section 4) tends to zero at late times if and only if 7 < | . 

1 From now on we will refer to this paper as WHU 

2 Prom a dynamical systems point of view, the late time behaviour signifies that the evolution is governed by the future 
attractor of the dynamical system. 
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In that paper a complete proof of the late time behaviour of dust models (in fact all models with | < 7 < |) 
was given, but the results for the radiation models were stated without proof. The goal of the present paper, 
which should be read in conjunction with WHU, is to remedy this deficiency. 

The late time behaviour of the Bianchi VIIo radiation models was first discussed by Doroshkevich et 
al (1973) (see equation (9) and equations (1.14) and (1.15) in Appendix I), and subsequently interpreted 
in a series of papers by Lukash (1974, 1975, 1976). Doroshkevich et al gave the asymptotic forms of the 
metric coefficients, but did not give a proof or any justification of their results. Reference is made to two 
preprints that we have been unable to locate (their references 19 and 20). Our approach, which makes use 
of the orthonormal frame formalism and expansion-normalized variables, provides a convenient framework 
for proving the results, and is also well suited for numerical simulations. In addition, it has the advantage 
that it leads directly to the asymptotic behaviour of the dimensionless anisotropy scalars (see Section 4) that 
describe the deviation of the models from a Friedmann-Lemaitre (FL) model. 

The outline of the paper is as follows. In Section || we present the evolution equations for Bianchi VIIo 
models with a radiation fluid. In Section |^ we establish the limits of the dimensionless gravitational field 
variables at late times. Section |i| is devoted to the question of isotropization and whether the Bianchi 
VIIo models are close to any FL model at late times. In Section || the detailed asymptotic forms of the 
gravitational field variables are found using centre manifold theory. In Section ^, the results of the previous 
Sections are used to determine the asymptotic forms of the metric components, and a comparison with the 
results of Doroshkevich et al (1973) is made. We conclude in Section (7) with a discussion of the implications 
of the results. The details of the proofs are given in three appendices. 



2 Evolution equations 

Since the evolution equations for Bianchi VIIo models with a radiation-fluid source can be found by specializ- 
ing the equations given in WHU, we will only briefly discuss the choice of variables. We use the orthonormal 
frame approach of Ellis & MacCallum (1969). The basic variables are {H, <r±,n±}, where H is the Hubble 
parameter. The variables a±, which parametrize the non-zero components of the shear tensor, describe the 
anisotropy of the expansion of the cosmological model, while the n± describe the spatial curvature of the 
orbits of the Bianchi VIIo symmetry group. The next step is to introduce expansion-normalized variables 
according to 

{S ± ,V ± } = l{ ( r ± ,n ± } , n = ^j, (3) 

(see Wainwright & Ellis (1997), page 112, for the motivation of this normalization). The density parameter 
f2 is given by 

= 1 - T,% - E 2 _ - , (4) 

(see equation (3.9) in WHU). The assumption of a non-negative fluid energy density yields O > 0. This, in 
conjunction with (|4|), implies that the variables £+,£_, and iV_ are bounded. The remaining gravitational 
field variable, N+, however, need not be. In fact, it has been shown in WHU that 

lim N + = 00 , (5) 

T — >00 

for all initial states with > and | < 7 < 2. The time-variable r is the so-called Hubble-time, which 
satisfies 

fr = ' (6) 

where t is the cosmological clock time. The behaviour of N + at late times motivates the introduction of a 
new set of variables {R,ip,M} according to 

= Rcosip , iV_ = Rsinip , M — , (7) 
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where R > 0. The case R — 0, which results in £_ = 7V_ = 0, describes the Bianchi VIIo locally rotationally 
symmetric (LRS) models. For these models, ip is irrelevant, as can be seen from (Q). Since Bianchi VIIo 
LRS models are equivalent to Bianchi type I LRS models, they are asymptotically self-similar and do not 
exhibit Weyl curvature dominance at late times. We will therefore not consider this particular class of models 
further. The evolution equations of the gravitational field are 

= -R 2 - ^+(1 - £+) + (1 + S+)i? 2 cos2V' , (8) 
R' = [(1 + £+)£+ + (R 2 - 1 - E+)cos2?A] R , (9) 
M' = - [(1 + S+) 2 + R 2 (cos 2^ + 3M sin 2^)] M , (10) 

f = i[2+(l + E + )Msin2^] , (11) 

where a prime denotes differentiation with respect to r. These equations correspond to equations (3.17)- 
(3.22) in WHU with the specific value 7 = |. Equation (||) now takes the form 

fl = 1 - Y? + - R 2 , (12) 

and it follows from (||) and (^|) that 

fl' = 2 + i? 2 cos 2tp) fl . (13) 
For future reference, we introduce an auxiliary variable Z according to 

which satisfies 

Z' = [-1 - (2 + S+)cos2?A-3i? 2 A/sin2V'] 2 , (15) 
as follows from (|),@, and (||). 

3 Limits at late times 

In this Section we find the limits of E + and R as r — * 00. The limit of M is 

lim M = , (16) 



which follows immediately from (|5[). Since ip' in ( |ll| ) is proportional to M 1 , the trigonometric functions 
111 (I)-© will oscillate increasingly rapidly as t — > +00 . To facilitate the study of the behaviour at late 
times, we therefore introduce new gravitational variables, fl+,R, and M, according to 

£+ = £+ - \M{\ + E+)i? 2 sin2V' , (17) 

_ (18) 

1 + \M(R 2 - 1 - S + )sin2V' ' 

M 

M = . (19) 

We likewise replace the auxiliary variables Z and f2 by the variables Z and f2, defined by 

^ = i t^- ; , (20) 

l + ±Af(l + £ + )sin2V> 

& = i— —o ■ (21) 

1 + ±Afi? 2 sin2?/> 
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These new variables are denned so as to ensure that the oscillatory terms in their evolution equations tend 
to zero as r — > +00 . This follows since the oscillatory terms in (p|)-(^0|) are bounded as r — > +00 . The 
evolution equations, which follow from (Iq)-(O), (|l3|), and (fTa), have the form 



S' + = -R 2 - £+ (l - £+) + MR 2 B t+ , (22) 

i?' = [E + + t\ + MB R ] R , (23) 

M' = -[(l + E + ) 2 +M%]M, (24) 

2' =[-1 + MB 2] 2 , (25) 

fi' = 2 [E 2 + Mi? 2 S n ] Q , (26) 



where the B's are bounded functions as r — > +00 . Equations (22)-(|26|) depend explicitly on -0 through the 
B's which are all multiplied by a factor of M. This fact, in conjunction with ( |16|) , implies that the oscillatory 
terms in (p2[)-(p6|) tend to zero as r — > +00, as desired. The structure of (^2j)-(|2q) now allows the limits 
of £+ and R as r — > +00 to be found, and subsequently the corresponding limits of E + and R. The main 
result is contained in the the following theorem. 

Theorem 1 

Any solution of the differential equations (p|)-(pH) with > that is not LRS, satisfies ( pH| ) and 

lim S+ = , lira -R = , 

r — >+oo r — >+oc 

with 

M „ 
hm — = . 

T^ + OC R 



Proof: The proof, which is based on equations (|22|)-(|26|), is given in Appendix A. 
■ 

In the next Section we will consider the physical implications of this theorem as regards isotropization 
of the cosmological model. 



4 Anisotropy of the models at late times 

To quantify the deviation of non-tilted Bianchi VIIo models from the flat FL model, we introduce a set of 
dimcnsionless anisotropy scalars, {E,iS, W}, all of which are zero for any FL model. The anisotropy in the 
expansion of a cosmological model is described by the shear scalar S, defined by 

_ —ah 

^ = ' (27) 

where o~ a b is the rate of shear tensor of the fluid congruence^. The anisotropy in the spatial geometry is 
described by the scalar S, defined by 

3 c 3 qab 

S 2 = » (28) 

where 3 S a b is the trace-free spatial Ricci tensor (see Wainwright & Ellis (1997), page 29). Finally, the 
intrinsic anisotropy of the gravitational field is described by the scalar W, defined by 

W 2 = E ^ Eab + H ^ Hab , ( 29 ) 

where E a b and H a b are the electric and magnetic part of the Weyl tensor with respect to the fluid congru- 
ence respectively (see Wainwright & Ellis (1997), page 19). The scalars S and W describe the dynamical 

3 It is this tensor that is parametrized by the variables <j± in (ph. 
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importance of the spatial curvature and the Weyl curvature compared to the overall expansion of the uni- 
verse respectively. The use of expansion-normalized anisotropy scalars dates back to the important paper 
on cosmological observations by Kristian & Sachs (1966) (see page 398), where they showed that various 
geometrical quantities, including the Weyl tensor can, in principle, be restricted by observations. For non- 
tilted Bianchi models, a zero shear tensor, — 0, implies that the Weyl tensor is zero, and can thus be 
viewed as characterizing the FL models. Restricting the expansion-normalized shear scalar E to be small, 
however, does not guarantee that the expansion-normalized Weyl tensor, as described by W, is small. For 
a model to be close to an FL model as t — > +oo , we therefore require that 

E < 1 , S < 1 , W < 1 , (30) 

as t — > +oo . The limits of the anisotropy scalars {E, S, W} asu +oo for non-tilted Bianchi VIIo models 
with a radiation fluid are given by the following theorem. 

Theorem 2 

For any non-tilted radiation-filled Bianchi VTio cosmological model that is not LRS, the density parameter 
satishes 



and the anisotropy scalars satisfy 



lim n = 1 

T — » + 00 



lim E = , lim W = +oo , 

r — >+oo r — >+oc 



and 

lim sup S — +00 , lim inf 5 = 0. 

Proof: This theorem is an immediate consequence of Theorem 1, once E, S, and W are expressed in terms 
of the variables E + , R, M, and i/j. In WHU it was shown that 

E 2 = E+ + R 2 cos 2 ip , (31) 

and 

W=™[l + 0(M)} , (32) 

(see WHU, equations (3.33) and (3.39)). Since 3 S a b is diagonal for Bianchi type VIIo models, it follows that 

S 2 = S 2 + + S 2 _ , (33) 
(see Wainwright & Ellis (1997), page 123), where 

s+ = Jh ^ 22 + 3 ^ 33 ) ' s - = W {3S22 3S33) • (34) 

Specializing equation (6.36) in Wainwright & Ellis (1997) to Bianchi VIIo models yields 

S+ = 2N-N+ S- = 2N 2 , (35) 
which, in conjunction with (0), implies that 

[! + £>(*•)]. (36) 
Equations (|l|), ©, and @ and Theorem 1 now give Theorem 2. 



The limits in Theorem 2 show that radiation filled Bianchi VIIo models isotropize with respect to the 
shear, but not with respect to the spatial curvature or the Weyl tensor. In view of (30), therefore, a general 
non-tilted Bianchi VIIo model with a radiation-fluid source is not close to any FL model at late times. 
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5 The asymptotic solution at late times 



In this Section we determine the asymptotic form of the gravitational field variables E_|_, R, M, and ip as 
t — > +00 . From these expressions, the asymptotic forms of the anisotropy scalars of the previous Section 
can subsequently be found. Theorem 1, in conjunction with (p4|), implies that 



M = O (e- bT ) 



(37) 



as t — > +00 , where the constant b can be written as b = 1 — 6. Here the constant S > can be chosen arbi- 
trarily small. It now follows from (|25| ) that there exists a constant Cz such that Z = Cz^T T [l + ( e ~ br )] • 
Hence, using d2fl), we find 



z = c z c~ T [l + O 



(38) 



as t — > +00 . The asymptotic behaviour of E + and i?, and hence of E + and R, is governed by (|22j) and 
( |23|). W e begin our analysis by considering the differential equations in M 2 obtained by omitting the terms 
~~ 1) that involve M. The resulting differential equations are 



in <\2 2h and 



# = E+(l + E+)i? , 



si 



-E+(l 



i? 2 



(39) 
(40) 



which we will refer to as the truncated differential equations. Their solutions will be denoted by E + and 



R. Since the proof in Appendix A that any solution of ( |22| ) and (|23|) with R > and > satisfies 
lim r _ 
that 



■+00 



E + = = lim r ^ +oc R also applies to the truncated differential equtions ( |39| ) and (|io|), it follows 



lim E+ = 



lim R = 

' — >-\-oo 



(41) 



for any solution with R > and il = 
for the truncated equations (|39|) and 



f 



E^ — i? 2 > 0. In other words, the origin is the future attractor 
A straightforward calculation shows that this equilibrium point 
: —1 and \fi , = 0. It will therefore be necessary to consider centre 

manifold theory^ when determining the asymptotic forms of E + and R as r — > +00 . The results are 
summarized in the following theorem. 



is non-hyperbolic, with eigenvalues Ag + 



Theorem 3 

For any solution of the truncated differential equations ( |3^ ) and ( fl(| ) with > 0, and R > 0, there exists a 
constant Cr such that 



R = 



/2t 
+ 2r 



1- —(lnr + C R )+0(T 2 
At 



l- — {\n T + C R -2) + 0(T 2 ) 



(42) 
(43) 



where T = — . 

T 

Proof: The proof is given in Appendix B. 



The next step is to relate the solutions of (E2J) and (g3|) to the asymptotic forms in Theorem 3. An important 
consequence of Theorem 3 is that E + and R exhibit a power-law decay to zero. This is the typical behaviour 
for non-hyperbolic equilibrium points. Since (|2^) and ( p3| ) differ from the truncated differential equations 
( |39| ) and ( |40| ) by the addition of exponentially decaying terms, it is plausible that the asymptotic decay rates 
of the solutions will not be affected. We therefore conjecture that any solution (R, E+) of ( p^ ) and ( p3| ) with 

4 For an introduction to centre manifold theory, see, for example, Carr (1981). 
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R > and 51 > will have the same asymptotic form as the solutions (R, E + ) of the truncated differential 
equations (39) and (ES) as given in Theorem 3, i.e., 



R = 



1 


1 - 


1 




47 


l 


1 - 


1 


~27 


2~7 



0(T 2 ) 



(44) 
(45) 



The dependence of the constant Cr on the initial conditions will be different than that of the solution (f42|)- 
( |43| ) . We have not been able to give a proof of the above conjecture in general, but in Appendix C we give 
a proof of a related but simpler result. Assuming ([5]) and (|45|), it follows from ([TJ]), and ( |37| ) that 



i? = 



The asymptotic form of M can be obtained from ( [46|) and 

1 



1 

V27 


1 

1 

4r 


(lnr- 




l 

~27 




(ln rH 


-Cr- 



M = C M Te~ T 



1 



4r 



2) + O (T 2 ) 
) using ([IJ) and 
(21nr + 2C^-l) + 0(T 2 ) 



The result is 



where C 



2Cz- Subsequently, by using (|11|), we find 
2e 



[l-^(21nr + 2CH-l) + 0(T 2 )] 



(46) 
(47) 

(48) 

(49) 



This equation determines tp up to an additive constant C^, although tp cannot be expressed in elementary 
terms. 

The Hubble variable H can be determined algebraically through the relation fl = jgr, since, for a 
radiation fluid, the quantity fil 4 is a constant. Here / is the length scale function, which is related to r 
according to 



I = Qe T , 



(50) 



where C; is a constant (see WHU, equation (1.3)). It follows, using (^), (|46|), ([47|), and ( p0[ ) that H is given 

by 



H = C H e~ 



1+ + 0(ter) 



(51) 



where C/f is a constant. Finally, by using (pl|), in conjunction with (^|), the relationship between the 
cosmological clock time t and the Hubble time r can be found. By choosing appropriate initial conditions 
for t, we obtain 



t = 



2C 



H 



(52) 



as r — * +oo . The freedom to shift the time- variable according to t — > r + C, where C is a constant, can 
be used to fix one of the constants, say C/ = 1. There thus remain four essential constants of integration, 
namely Cr,C^,,Cm, and Ch- This confirms that (f4(f)-(|5"2|) describe the asymptotic form of the general 
Bianchi VIIo solution with a radiation fluid source. 

We conclude this Section by giving the asymptotic form of the various scalars of physical interest in terms 
of the clock time t. The asymptotic form of the anisotropy scalars S 2 , W 2 and S 2 follow immediately from 
(|3l|), (||^), (^)-(|49|) and (|52]) . The asymptotic form of the Hubble variable H and the radiation density /i 
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are obtained from (|J), ([l2|), (|46|), (f47j), ( pT{ ) and (|52]). For convenience we introduce new constants A and a 
according to 



A = 2C H , a = — . (53) 



Note that A has dimension of (time) , while a is dimensionless. The results are 



^( 1+ 21nTII)) < (54) 

' (55) 

2 _ cos 2 y> 1 f56l 

E ~ ln(At) + [l n (^)] 2 ' (56) 

W 2 « , (57) 
[ln(At)f 

S 2 w W 2 sin 2 ^ , (58) 



where 



di ~ 21n(Ai) 



(59) 



It is understood that these asymptotic forms are valid for values of A and t such that \n(At) 3> 1. As 
mentioned in Section |l|, these asymptotic forms could describe the dynamics of the universe in the final 
stages of the radiation-dominated epoch, which ends at a time t = t oq , when the matter density first equals 
the radiation density. In order to describe this situation, the constant A must satisfy 

ln(At cq ) » 1 . (60) 

It follows from (|56|), in conjunction with (p9|), that the constant a determines how rapidly S oscillates. It 
also determines the magnitude of the Weyl scalar VV at t = t cq (see (|57|)). We shall discuss the significance 
of the above results in Section 7. 



6 Relation with the metric approach 



The late time behaviour of Bianchi VIIo models with a radiation-fluid source has been considered by 
Doroshkevich et al (1973) using metric variables. In our approach, the asymptotic behaviour of the metric 
components can be found directly from the behaviour of the orthonormal frame variables in (^6|)-(|5l|), and 
we now do this for the purpose of comparison. 

For a Bianchi type VIIo model, a set of group-invariant and time-independent 1-forms, u) a , a. = 1,2, 3 
satisfying 



du} 1 = , duj 2 =w 3 Au' , duj 3 =w 1 Au 2 
can be introduced for which the line element can be written 

ds 2 = -dt 2 + g a pLJ a LJ P . 



(61) 



(62) 



Here t is (as previously) the cosmological clock time. For a non-tilted Bianchi VIIo model, g a p can be 
diagonalizcd. Using the variables of Misner (1969) to parametrize g a p, the line element (62) can be written 
as 



ds 2 



-dt 2 + I 2 



-4/3" 



(63) 
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The relationships between the metric variables (3^ and the orthonormal frame variables for Bianchi class A 
models are given in chapter 9 of Wainwright & Ellis (1997). Specializing these relationships to Bianchi VIIo 
models yields 



HI 

Isolating (3 + and f3~ leads to the following expressions 



HI 



l 2 H 2 

tanh(2\/3/r) = V3RM sin?/> , e 4/3 = - 3i? 2 M 2 sin 2 V) 

(see Appendix C in WHU). Using ©,©,©, and ©, it follows that 

1 



/3 + 
/3" 



■i^--lnr-^(Cfl-l), 



2V2 



l + ilnr+i^-l) 



sin^ . 



(64) 
(65) 

(66) 
(67) 



To facilitate a comparison with the asymptotic expressions of Doroshkevich et al (1973), we note that their 
variables 7, \i and (A1A2) 1 / 2 are related to the variables used in this paper by 



7 = I 6 , 11 = 4V3/T , (AiA 2 ) 1/2 = I 
and that their time variable T (the so-called BKL time) is related to r by 



dT 



Hr 



T«-(C fl C,)- 1 e- r 1 



4r 



(68) 



(69) 



We note that T approaches zero from below as r 
asymptotic forms: 



« £>i 



J^2 



13- ^ D 2 e{T)- 1 / 2 Tsmi>(T) , 



00 . Doroshkevich e£ aZ (1973) give the following 
dV> D 3 9(T) 



dT 



rp2 



(70) 



see Appendix I, equation (1.14). Note that we have relabeled their constants. The function 0(T) in Doroshke- 
vich et al (1973) is given by 



6(T) 



In v 

^7 



where v = \n(D 5 /T) (see Appendix I, equation (1.15) ). Using ( |69| ) in ( |7l"| ) yields 



.11-^- 



lnr 
~2t 



(71) 



(72) 



where the constant is a combination of the constants Ch , Ci, D4 and D§. Using ( |69] ) and (72) in ( |70| ) 
yields @ and @ if we identify C fl = 2L>„, £>i = C H 2 C*f 4 , and L> 2 = -V6/-D1. From @ it also follows 
that the evolution equation for tp in (70) is (asymptotically) equivalent to (EjJ) if we identify D3 = 2/D\. 



7 Discussion 



In this paper we have given a complete analysis of the behaviour at late times of non-tilted Bianchi VIIo 
cosmologies filled with radiation. In particular, we have established rigorously the late time behaviour of the 
dimcnsionless anisotropy scalars £,W, and iS, that determine the extent to which the models deviate from 
an FL model (Theorem 2). In addition, we have derived the asymptotic form as r — > +00 of the general 
solution, although in this regard we have had to rely on a heuristic argument at one stage of the analysis. 

Our results show that the Bianchi VIIo radiation models differ from the dust models in two significant 
ways: 
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1. for dust models, the dimensionless Weyl scalar W has a finite non-zero limit (see Theorem 2.4 in 
WHU), while for radiation-filled models, W increases without bound, 

2. for dust models, the dimensionless shear scalar £ tends to zero at an exponential rate in terms of the 
Hubble time r, while for radiation models £ decays at a power law rate in terms of r. 

The results also shed new light on the dynamics of the radiation models at late times. These models 
isotropize as regards the shear (i.e. £ — > 0), but not as regards the Weyl curvature and the anisotropic 
spatial curvature (i.e. W 0,6> 0). Equations (p4])-(p8|) give the asymptotic dependence of the various 
scalars of physical interest on the clock time t. We note that the leading order terms of H and /i in ( p4| ) 
and (]5^) give the flat FL radiation solution, and that, in addition, the shear scalar £ in (|5^) tends to zero. 
The behaviour of W and S in ([57]) and (|58|), however, implies that the solution is not close to the flat FL 
model at late times. Moreover, since W diverges at late times, the models exhibit the phenomenon of Weyl 
curvature dominance, i.e., the Weyl curvature plays a dominant role in determining the dynamics. These 
considerations show that solving the linearized Einstein field equations can lead to incorrect conclusions 
regarding the dynamics of the full equations. It was pointed out in WHU (see pages 2588-89) that some 
aspects of the late time evolution of dust models of Bianchi VIIo arc correctly described by linearizing the 
Bianchi VIIo evolution equations about the flat FL model. On the other hand, for radiation models, the 
linearized solutions totally fail to describe the asymptotic behaviour of the shear, the Weyl curvature, and the 
anisotropic spatial curvature at late timesf^ implying that this behaviour is a fully non-linear phenomenon. 

The rate at which the shear scalar £ tends to zero requires further comment. Equation ( |56| ) shows that 
the asymptotic form of £ 2 is the sum of two terms. One of them is oscillatory and decays as (in At)^ 1 , 
while the other is non-oscillatory and decays as (\nAt)~ 2 . It appears that the dominant oscillatory term 
has been inadvertently omitted in the literature giving the incorrect result that £ ~ \ n (At) as * 00 ■ 
For example, Doroshkevich et al (1973) state that "the anisotropy of the deformation tensor, decreases 
very slowly, proportionally to (ln£) -1 " (see page 741), but do not specify what normalization was used to 
calculate the "anisotropy of the deformation tensor". In a subsequent article, however, Novikov (1974), 
states the conclusion more explicitly, namely that the components of the shear tensor, when normalized with 
the Hubble scalar H, decay as [In^/tp)]^ 1 , (see equation (2) on page 275, and the definition^ of AHi/H on 
page 274). This incorrect result was subsequently quoted by Barrow (1976) (see the equation for a/6 in the 
case p — , on page 364). 

As mentioned in Section 1, the late time behaviour if Bianchi VIIo radiation models could describe the 
dynamics of the universe in the final stages of the radiation dominated epoch. Our results show that on 
transition from the radiation-dominated epoch to the matter-dominated epoch at t = t cq , the shear scalar 
would be small but the Weyl curvature scalar W would be large. These values of £ and W would then act 
as initial values, so to speak, for the matter-dominated epoch. Since £ — > and W — > Woo, where Woo is 
a positive constant, as t — > oo for a dust model, the shear scalar would remain small, but the dynamical 
significance of W would persist. It would be of interest to investigate density perturbations in such a model, 
since they could well behave differently than in a perturbed FL model. 
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Appendix A: Proof of Theorem |I] 

In this Appendix we give the proof of Theorem [lj The proof is based on the evolution equation for 0, and 
the fact that f2 is bounded. We begin by integrating equation (26) to obtain 

"0 Jr a 

where Clo = O(tq) and the initial time tq is chosen such that M < 1 for all r > tq. Equation ( |25| ) implies 
that 



\ In -2- = I t 2 + ds + I MR 2 BqcIs , (73) 



Z = oU- 1+5 ^ T ) , (74) 



where <5 can be chosen arbitrarily smallQ We now require the fact that Q is bounded, which follows from 
(f|), @, and (|l|). Recalling the definitions (0) and @, equation @ implies 

Mi? 2 = O (e(- 1+ ^ T ) , (75) 

and hence 

f MR 2 B n ds = C + oU- 1+ ^ T ] , (76) 

where C is a constant depending on the initial conditions at r = To. Using ( |76| ) and the fact that fl is 
positive and bounded, we conclude from ( f73| ) that 7(r) = S^ds is bounded above for all r > t . Since 
the integrand is non-negative, the function Z(r) is monotonically decreasing, which implies that lim r _,_|_ 00 I(r) 
exists and is finite. It follows from (|73|), ( |76| ) and (|2l]) that lim T ^ +00 f2, and hence lhxii-^+oo Q exists, i.e. 
there is a constant L, satisfying < L < 1, such that 



lim O = i . (77) 

T — >00 

In order to proceed, we need the following standard result 
Lemma 1 

Let f(r) be a non-negative real-valued and Lipschitz continuous function on the interval tq < r < +oo. If 

f(r)dr 

TO 

is finite then lim r _> +00 / (t) = 0. 

Since £+ and R are bounded, it follows from ( p^ ) that (£+)' = 2£+£^ is uniformly bounded for r > tq. 
This implies that £+ is Lipschitz continuous on the interval r > to (see ???). We may now apply Lemma [l], 
with / = E+, to conclude that 

lim E+ = , (78) 

T — > + 00 



and hence, from ( |l6| ) and ( p7| ) 



lim £+ = . (79) 



From m,(m,(Wn, and (79), it follows that 



lim R 2 = lim R 2 = 1- L . (80) 

T^ + OC r— > + OC 



7 Since M tends to zero and is bounded, one can choose tq such that |AfSi| < S for r > tq. 



n 



Equation (p2|), in conjunction with J79|), (jS0|), and (|l6|), now implies that lim T _^ +oc exists, and that 

lim £+=£-!. (81) 



T — > + OC 



This result contradicts (|78|) unless L = 1. Using ( |80| ) we obtain 

lim i? = . (82) 



r — >+oo 



To find the limit of the quotient ^ , we note that (|2^) and ( pi[ ) can be combined to give 



-g) = (-1 - 3S+ - + MS*) I — ) , (83) 



where £?* is a bounded function as r — > +oo . It now follows from (f7q) and (161) that 



and hence that 



M 

lim --=0, (84) 

T-> + 0O R 



M 

lim -=- = 0. (85) 

T^ + OO K 



Appendix B: Proof of theorem 3 

In this Appendix we give the proof of Theorem 3, using centre manifold theory. In developing this theory, 
Carr (1981) (see his equation (2.3.1)), considers a system of differential equations of the form 

x' = Ax + f(x,y) , (86) 
y' = By + g(x,y) , (87) 

where x £ R™, y £ M. m , A is a constant n x n matrix, i? is a constant mx m matrix, and / : R™ x R m — > R™ 
and g : 1" x R m — > R m are C 2 functions that are zero and have zero derivatives at the equilibrium point 
(x, y) = (0, 0). It is assumed that the eigenvalues of A have zero real parts and that the eigenvalues of B have 
negative real parts. The assumption on A means that the equilibrium point (x, y) — (0, 0) is non- hyperbolic, 
and hence that the linearization of (^) and (|87j) at this equilibrium point does not describe the behaviour 
of solutions near the equilibrium point. 



By theorem 1 in Carr (1981), the differential equations (|86|)-(87) has a local center manifold y = h(x), 
where h £ C 2 , and the flow on this center manifold is governed by the reduced differential equation 

u' = f(u,h(u)) , (88) 

(see equation (2.4.1) in Carr (1981)). By theorem 2 in Carr (1981) if the equilibrium point u = of ( |88| ) is 
asymptotically stable then the equilibrium point (x, y) = (0, 0) of (|8^)-(|87|) is also asymptotically stable. In 
addition, if (x(t),y(t)) is a solution of (j36|)-(lS7|) with (x(0), j/(0)) sufficiently close to the origin, then there 
exists a solution u{t) of ( |88| ) such that[] 

x(t) = u{t) + O (e~ 7 *) , (89) 
y(t) = h(u(t)) + O (e-^*) , (90) 



8 It is possible that the solution u(t) is the zero solution, in which case x(t) and y(t) decay exponentially to zero. This special 
case will occur only if_thc initial point (x(0),y(0)) lies on the stable manifold of the equilibrium point (x,y) = (0,0) of the 
differential equations (p6|)-(|87|). 
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where 7 > is a constant depending only on the matrix B (see Carr (1981), equation (2.4.5)). The implication 
of this result is that if the equilibrium point u — of flSq) is asymptotically stable, one can determine the 
asymptotic behaviour of the solutions of (pq)-(|87|) by finding the asymptotic behaviour of solutions of (pS|). 
In order to perform this analysis one has to approximate the center manifold y — h(x), as follows. For a 
function <fi : l n — > t m that is C 1 near the origin, define an operator M by 

(M0) (x) = <f>'{x) [Ax + f(x, - + g(x, <f>(x))] . (91) 

By theorem 3 in Carr (1981), if 0(0) = 0, 0'(O) = and [M<j>) (x) = O (\x\ q ) as x -> 0, where g > 1, then 

\h(x) - (t>(x)\ = O (\x\ q ) , (92) 

as a; — ► 0. In other words, the function <p(x) approximates the center manifold y = h(x) with an accuracy 
determined by q. 

The differential equation ©-© is of the form @-@ with n = 1, m = 1, 



(x,y) = (R, E+) , 



(93) 



A = , B = 1 , 



(94) 



= a:y(l +y) , y) = -x 2 + y 3 



(95) 



In addition, we have shown that the equilibrium point (i?, £+) is asymptotically stable, so we can apply 
the above theory. Using (|9l]) and (|92|), in conjunction with (^3|)-(^), we find that the center manifold is 
approximated by 

£+ = -R 2 -2R 4 + (r 6 ^) . (96) 
Substituting this in (|39| ) gives the reduced differential equation ( |88| ) on the center manifold in the form 

R' = -R 3 -R 5 + (r 7 ) . (97) 



This differential equation is precisely equation (3.1.9) with a = —2, considered by Carr (1981), who shows 
that the asymptotic form of the solution is 



R = 



1 



>2t 



l-i-(lnr + C fi ) + (r- 1 ) 



(98) 



as r — > +00, where Cr is a constant that depends on the initial conditions. The error bound o(r x ) can be 
improved to O ^(~^-) 2 ^ ; as hi Theorem 3, by increasing the accuracy of the center manifold approximation 

£ + = -R 2 -2R i -llR 6 + o(R 8 ^ . (99) 
Finally, substituting j9q) in (p6h gives the asymptotic form of S+, as stated in Theorem 3. 



Appendix C 

In this Appendix we prove that any solution to the perturbed differential equation 



-y 3 + ee-~> T 



(100) 
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where e > and 7 > are constants, has the same asymptotic form as r — * +00 as the solutions to the 
unperturbed differential equation 



Y > = _y3 



(101) 



Equation (101) is a simplified form of the differential equation (|97| ) that describes the flow on the center 
manifold in Appendix B, and the term ee~ 7r mimics the perturbation introduced by the exponentially 
decaying function M. 

Proposition 1 

For any solution y(r) of \10(\) there exists a solution Y(t) of such that 



as t — ► +00 



|y(r)-F(r)|=0(e-^) , 



Proof: Let y(r) be the solution of (10C) that satisfies the initial condition y(0) = y$ where yo > 0. The 
differential equation (100) then implies that y(t) remains positive and is bounded above. This, in turn, 
implies that the solution exists for all r > 0. 



We first need to derive an upper bound for ^pj- The differential equation (10C) implies that 



which, when integrated from to r, gives 



We now define 



V 



1 „ 1 

-T^ < 2T + — 



(102) 



(103) 



w(t) 



The differential equation (100) implies that 



Mr) 



ee 



Integrating fllOq) from to r gives 



w(t) — w(0) = —e 



y( s ) : 



:d,S 



Since is bounded (see (103)), it follows that lim T ^ +00 w(T) exists and equals c, where 



c = w(0) — e 



y( s ) 3 



ds 



Subtracting ( |107 ) from (^06j) and using ( |103| ) yields 



w(t) — c = e 



r y(s) 3 



<ls : . / i 2.s - -5 I <ls 



1 



.3/2 



Thus, 



;(r) = c + O (t 3/ V 7T ) 



(104) 



(105) 



(106) 



(107) 



(108) 



(109) 
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as r — + +00 . It now follows from (104) that 



y(r) = , = , (no) 

v^Jt + c + O^/V^) 



which can be rearranged to give 



y(r)= 1 +0(e~^) , (111) 

V2y/T + C 



as r — * +00 . Since the first term is a solution of (101), the proof is complete. 
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